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Using molecular dynamics simulations we examine the effective interactions between two like- 
charged rods as a function of angle and separation. In particular, we determine how the competing 
electrostatic repulsions and multivalent-ion-induced attractions depend upon concentrations of sim- 
ple and multivalent salt. We find that with increasing multivalent salt the stable configuration of 
two rods evolves from isolated rods to aggregated perpendicular rods to aggregated parallel rods; 
at sufficiently high concentration, additional multivalent salt reduces the attraction. Monovalent 
salt enhances the attraction near the onset of aggregation, and reduces it at higher concentration of 
multivalent salt. 
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Multivalent-ion-induced aggregation of stiff polyelec- 
trolytes has been studied extensively in recent years for 
at least two important reasons: it arises from correlations 
not included at the mean-field (Poisson-Boltzmann) level 
0, H H II U 0, 0, H and it lies at the heart of bi- 
ological phenomena such as cytoskeleton re-organization 
[Hi [H 13 and DNA packaging [H 03. Most stud- 
ies have focused exclusively on the interaction between 
charged rods that are parallel. A few theoretical calcula- 
tions have allowed for non-parallel orientations 0,0], 
but only at a single multivalent salt concentration. Re- 
cent work jl 7L however, suggests that the evolu- 
tion of angle-dependent attractions with changing multi- 
valent salt concentration is crucial to phase behavior. For 
example, at intermediate multivalent salt concentrations 
F-actin solutions can form lamellar phases of stacked 
rafts where each raft consists of two layers of mutually 
perpendicular actin filaments, but at higher multivalent 
salt concentrations they form bundles of nearly parallel 
filaments 01- Similar physics may apply to transitions 
from networks to bundles of F-actin in the cytoskeleton 

m. 

In this study we use explicit-ion, continuum-dielectric 
molecular dynamics (MD) simulations to examine how 
the concentrations of monovalent and multivalent salts 
affect the angle- and distance-dependent effective poten- 
tial between two charged rods. A threshold concentra- 
tion of multivalent salt is needed for the two rods to at- 
tract. We find that the preferred configuration is perpen- 
dicular just above the threshold, and parallel at higher 
concentrations of multivalent salt, in agreement with the 
experiments on F-actin solutions cited above 01- f ur " 
thermore, monovalent salt can lower the concentration 
of multivalent salt needed for attraction, suggesting that 
monovalent salt may induce aggregation. Finally, we 
find that the rods can still attract even when they are 
overcharged Q (i.e. the sign of their effective charge 
is reversed by condensed counterions) at high multiva- 
lent salt concentrations, and that the attraction weakens 



with increasing overcharging. However, we never observe 
strong overcharging because above a certain concentra- 
tion of multivalent salt the added multivalent ions si mply 
form complexes with monovalent co-ions in solution [2(1 ] . 

In our simulations each rod is composed of 64 spheri- 
cal monomers, each carrying a charge of -1 in units of the 
electronic charge e, separated at fixed intervals. The two 
rods are perfectly rigid and fixed at a specified center-to- 
center separation, R, and angle 7. In all cases, one rod 
lies parallel to the face diagonal of the enclosing periodic 
box I2J, while the other is rotated away from a paral- 
lel configuration at an angle 7 about the axis connecting 
the centers of the rods. In addition, we introduce mobile 
multivalent ions of charge +3, and mobile monovalent 
ions of charges +1 and -1. The system is always electro- 
statically neutral, with 128 ions of charge +1 to balance 
the charge on the two rods, one ion of charge -1 for every 
additional ion of charge +1 (monovalent or 1:1 salt), and 
3 ions of charge -1 for every ion of charge +3 (trivalent 
or 3:1 salt). 

We reference salt concentrations to the total charge on 
the two rods (128 electronic charges). For example, a 
3:1 salt concentration of 03a = 1 means that the total 
charge due to +3 ions is equal to the total charge on the 
two rods. Similarly, a 1:1 salt concentration of cia = 1 
means that the total charge due to +1 ions from the 
monovalent salt is equal to the total charge of the two 
rods. 

We include two types of pair interactions between par- 
ticles. First, we use the truncated Lennard-Jones poten- 
tial to allow for short-range repulsions. This introduces 
the energy scale e and the particle size a. Second, we 
include the Coulomb interaction, Z\Zij eryi, where Zi is 
the charge on particle i and e is the dielectric constant. 
To handle the long-range Coulomb interaction in our sys- 
tem with periodic boundaries we use the Particle Mesh 
Ewald (PME) method [12]. 

Our simulations are carried out in the canonical (NVT) 
ensemble with the temperature fixed at k-sT = 1.2e us- 
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ing the Langevin thermostat [23j. The monomer num- 
ber density, 10~ 4 cr~ 3 , corresponds to a box volume of 
(109er) 3 . The dielectric constant is chosen such that the 
Bjerrum length (the distance at which the electrostatic 
interaction between two electronic charges is equal to the 
thermal energy) is Ib = l/ek^T = 3.2ct. The separation 
between charges on the rod (i.e., the monomer separa- 
tion) is I = l.lcr and the dimensionless Oosawa-Manning 
ratio Q, Hi| is Ib/1 = 2.9, well above the threshold for 
counterion condensation. 

The characteristic time scale for the simulation is 
t = ^/mer 2 /e, where m is the particle mass. We use 
the leapfrog- Verlet integration scheme with a time step 
of O.OIt. Each simulation run is equilibrated (as mea- 
sured by the leveling-off of the energy) for at least 10 3 r 
before we collect data. The force on each monomer is 
then averaged over 40r intervals until we obtain 250- 
450 average-force data points. More extended runs were 
performed for a few systems to check that our results 
are not affected by the choice of equilibration time, time 
step or data collection interval. For each simulation run 
we calculate both the average normal force per monomer 
between the two rods and the average torque on the rods 
about the center-to-center axis. Error bars (indicated by 
the size of the points in our figures) correspond to the sta- 
tistical error associated with the average forces at each 
40r interval. 

To obtain the effective interaction potential, or re- 
versible work, as a function of the separation R between 
the two rods, we should integrate the normal force with 
distance as the rods are brought from infinite separa- 
tion to R at fixed angle 7. This is not possible within 
our periodic-boundary-condition simulation. Instead, we 
calculate the reversible work, AW(R), for bringing the 
two rods from a fixed reference separation of 8cr to R 
at fixed 7. Similarly, the effective potential as a function 
of angle, AW^), is calculated by integrating the torque 
from 90° to 7 at fixed R. 

We first ask how the effective interaction between two 
charged rods depends on the concentration of multivalent 
salt, in the absence of monovalent salt. Fig. ^a) shows 
a plot of the effective potential per monomer, AW, as 
a function of separation, R, when the two rods are par- 
allel (7 = 0°). In the absence of 3:1 salt AW(R) is 
positive for all R, implying that the rods repel. When 
enough 3:1 salt is added (c 3: i » 0.3), the effective poten- 
tial develops a global minimum at small R; the rods now 
attract each other. This is the threshold 3:1 salt con- 
centration for aggregation; note that this corresponds to 
approximately only one-third of the charge on the rods 
neutralized by trivalent ions. Beyond this threshold the 
attraction increases with multivalent salt concentration 
until C3 : i = 1; this is where the charge on the rods is 
completely neutralized by the trivalent counterions. Be- 
yond C3 : i = 1, the attraction decreases slightly (short- 
dashed curve; C3.1 = 12). Thus, the magnitude of the 



attraction is non-monotonic with 3:1 salt |2jj]. Similar 
non-monotonicity is observed in experiments with mixed 
salts on DNA solutions (13] ■ Though it is imperceptible 
in the figure, the effective potential beyond 4.5er is pos- 
itive for both C3:i = 1 and C3.1 = 12, which means that 
the rods still repel one another at large distance. 

The angle dependence of the effective rod-rod inter- 
action is shown in Fig. db) for R = 2.1a, near the at- 
tractive minimum. In the absence of 3:1 salt the rods 
prefer to be perpendicular since the rods repel one an- 
other. Just above the threshold for aggregation (solid 
curve), the global minimum is at 90°, implying that the 
preferred aggregated configuration is a cross. As the con- 
centration of multivalent salt increases further the min- 
imum at 7 = 0° deepens, and at c 3: i = 1 the rods now 
prefer to aggregate in a parallel configuration. When still 
more multivalent salt is added the minimum at small an- 
gle decreases; this shows again that the effective potential 
depends non-monotonically on the concentration of mul- 
tivalent salt. 

In many-rod systems under conditions for which a pair 
of rods prefers to aggregate perpendicularly, we would 
expect to find networks, rafts or other structures where 
rods cross each other at large angles. On the other hand, 
when a pair prefers to aggregate in a parallel configura- 
tion, we expect to find bundles or networks of bundles. 
Note that for a pair of rods, the preferred angle at which 
they aggregate does not vary continuously with increas- 
ing multivalent salt as it would at zero temperature, but 
rather jumps from 90° to 0°. This entropic effect suggests 
that transitions from network/raft structures to bundles 
may be first-order (l^ |. 

Fig.|2shows the effects of added monovalent salt on the 
effective interaction at a multivalent salt concentration of 
C3 : i = 0.3, just above the threshold for aggregation. In 
Fig. |2I a) the attraction is seen to grow stronger with in- 
creasing 1:1 salt. Furthermore, we find (not shown) that 
if the 3:1 salt concentration is below threshold, then 1:1 
salt can actually drive the pair to aggregate. Fig. Efb) 
shows that the metastable minimum at 0° increases in 
depth and becomes the global minimum with increasing 
1:1 salt. This result suggests that adding 1:1 salt to a 
solution of charged rods in a network or raft phase can 
drive the solution into the bundle phase. To our knowl- 
edge, an experiment has not yet been designed to test 
these two effects of added 1:1 salt. 

To see why monovalent salt enhances the effective at- 
traction, we measure the neutralization fraction, f(r c ), 
defined as the time average of the sum of all mobile 
charges, positive and negative, within radius r c from ei- 
ther one of the rods, divided by the total charge of a rod 
[ifij . The sketches in Fig. depict a top view of the vol- 
ume that encompasses the mobile charges included in the 
calculation of f{r c ) for a given r c and separation R. For 
r c < R/2 (to the left of the vertical dotted line), there 
are two separate enclosing volumes, as sketched, while for 
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FIG. 1: (a) The effective potential, AW(R), between two 
parallel charged rods for different concentrations of multiva- 
lent salt. The size of the solid circles at C3.1 =0 corresponds 
to the error bar for all points on all curves, (b) The effective 
potential, AW, as a function of angle, 7, for two rods sepa- 
rated by R = 2.1cr, for different concentrations of multivalent 
salt. 



FIG. 2: The effective potential between two charged rods, 
AW as a function of (a) separation R for parallel rods (7 = 
0°) and (b) angle 7 for closely-separated rods (R = 2.1<r), 
for different monovalent salt concentrations. The multivalent 
salt concentration is fixed at C3 : i = 0.3, just above the thresh- 
old for aggregation. The size of the solid circles at C\-\ = 
indicates the error bar for all points on all curves. 



r c > R/2 (right of the vertical line), the two enclosing 
volumes merge into one. Fig. shows /(r c ) for different 
concentrations of 1:1 and 3:1 salts. The solid (dashed) 
curves correspond to cases without (with) 1:1 salt. Near 
the 3:1 salt threshold for aggregation, c^-i = 0.3, the 
solid and dashed black curves show that f(r c ) increases 
with added 1:1 salt. We find that the effective charge 
on the rods is reduced by nearby counterions from the 
1:1 salt while the concentration of multivalent counteri- 
ons near the rods is nearly unchanged j^j. As a result, 
the repulsive contribution to the effective interaction be- 
tween rods is reduced while the attractive contribution 
is unaffected, giving rise to the increase in the net effec- 
tive attraction shown in Fig. (Note that for the black 
dashed curve, the screening length-about 6<r-is signifi- 
cantly larger than the inter-rod separation of R — 3.6c 
and hence not relevant here.) 

At high concentrations of multivalent salt (c^.i > 1), 
we find that monovalent salt has the opposite effect-it ac- 



tually reduces the effective attraction. This is consistent 
with experiments on DNA solutions [22]]. This effect also 
can be understood by looking at f(r c ). The gray solid 
and dashed curves in Fig. show that /(r c ) decreases 
with added 1:1 salt, consistent with earlier predictions 
[icj . Furthermore, we find that fewer multivalent ions 
contribute to f{r c ) upon the addition of 1:1 salt because 
the multivalent ions stay in solution in the form of com- 
plexes with monovalent co-ions. We conclude that co-ions 
lure multivalent ions away from the two rods, causing the 
attraction to decrease. 

The top curve (dotted) in Fig. shows that overcharg- 
ing occurs (f(r c ) > 1) at sufficiently high concentrations 
of multivalent salt. Note that/(r c ) crosses 1 at r c rs 1.3<t. 
Since R > 2r c at this point (the situation sketched to 
the left of the vertical dotted line) , the total charge en- 
closed in each cylinder shown is now positive and the rods 
are overcharged jjj. This overcharging saturates with 3:1 
salt; we find, for example, that when the multivalent salt 
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FIG. 3: Neutralization fraction / for 2 parallel rods at 
R = 3.6cr as a function of the radius of the encompassing 
volume r c for several different concentrations of multivalent 
and monovalent salts. The vertical dotted line corresponds to 
r c — R/2- to its left, we sum the mobile charges enclosed in 
two separate volumes, as sketched in top view, while to the 
right, we sum the mobile charges enclosed in a single volume 
with a figure-eight cross-section. 



concentration is tripled from C3.1 — 4 to C3.1 = 12 there 
is only a few percent increase in f(r c ). Beyond a cer- 
tain concentration of 3:1 salt, additional multivalent ions 
stay in solution in the form of complexes with oppositely- 
charged monovalent ions. As a result, we never observe 
strong overcharging. 

Fig. H shows that overcharging weakens the effective 
attraction. At first glance this is not surprising be- 
cause overcharging should increase the Coulomb repul- 
sion 0. However, when we increase the 3:1 salt concen- 
tration above C3.1 = 1.0 we find that the normal force 
between rods (which we integrate to obtain the effec- 
tive interaction) remains negligibly small at large sepa- 
rations, whereas it becomes significantly less negative at 
small separations. This implies that overcharging does 
not affect the contribution of the longer-ranged Coulomb 
repulsion but appreciably weakens the contribution of 
counterion-mediated attractions to the effective interac- 
tion. Fig.Elsuggests why the repulsion is not significantly 
affected by overcharging: f(r c ) is approximately as far 
above unity for C3.1 =4 as it is below unity for C3.1 = 1. 
Thus, the magnitude of the total charge enclosed within 
r c is the same in the two cases although the sign of the 
charge has flipped. 

How then does overcharging weaken the counterion- 
mediated attraction? As the concentration of multivalent 
salt increases, the amount of condensed charge increases 
slightly by accumulating primarily on the far sides of the 
rods. We find that this reorganization of charge (from 



"bonding" to "antibonding" regions) leads to a decrease 
of the attraction. 
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